Spontaneous potential (SP) logging remains a widely used method in well geophysics. However, its interpretation is often limited by simplified physical models and correction charts that do not fully account for the processes governing SP generation, particularly in shaly and heterogeneous formations. In this study, we develop a finite element-based algorithm for modeling SP responses in complex near-wellbore environments, with the aim of providing a more physically consistent framework for interpretation. The proposed algorithm is based on the numerical solution of the Poisson equation with electrochemical source terms, incorporating the cation transport number to describe diffusionadsorption processes and allowing for smooth variations in formation resistivity, fluid properties, and shale content. The numerical implementation is validated against published analytical solutions, correction charts, and previous numerical studies, showing good agreement in both the shape and amplitude of modeled SP responses across a range of geological scenarios, including thin beds and invasion zones. Application to real data from a Southeast Asia gas field demonstrates that the approach provides reliable estimates of formation water salinity and the cation transport number, with results consistent with independent estimates. The proposed method offers a flexible tool for SP response modeling and may complement existing interpretation techniques, particularly when working with heterogeneous formations and limited legacy datasets.